#### setting environment ####
require(lmtest)
require(sandwich)
require(clubSandwich)
require(effsize)
require(ltm)
require(dichromat)

#### data preparation ####
data <- read.csv("replication_data.csv")

## demographics and other characteristics
# respondent's gender
data$R.gender <- data$sex - 1
# respondent's age
data$age <- data$age + 17
# college of technology, junior college, or vocational school
data$middle.edu <- ifelse(data$education > 2 & data$education < 6, 1, 0)
# university or graduate school
data$high.edu <- ifelse(data$education > 5, 1, 0)
# political knowledge
data$knowledge <- (data$quiz1 == 1) + (data$quiz2 == 1) + (data$quiz3 == 4) + (data$quiz4 == 2)

# recoding missing data in the questions on political attitudes
data$nuclear.power[data$nuclear.power == 9] <- NA
data$defense[data$defense == 9] <- NA
data$constitution[data$constitution == 9] <- NA
data$public.security[data$public.security == 9] <- NA
data$separate.surnames[data$separate.surnames == 9] <- NA
data$preemptive.attack[data$preemptive.attack == 9] <- NA
data$same.sex.marriage[data$same.sex.marriage == 9] <- NA
data$income.tax[data$income.tax == 9] <- NA

## outcomes
data$DPRK.outcome <- ifelse(data$order == 0, data$experiment1, data$experiment2)
data$welfare.outcome <- ifelse(data$order == 0, data$experiment2, data$experiment1)

## experimental conditions
# female condition in the DPRK question
data$DPRK.female <- ifelse(data$DPRK.condition > 4, 1, 0)
# liberal condition in the DPRK question
data$DPRK.liberal <- ifelse(data$DPRK.condition == 3 | data$DPRK.condition == 4 | data$DPRK.condition == 7 | data$DPRK.condition == 8, 1, 0)
# aggressuve condition in the DPRK question
data$DPRK.aggressive <- ifelse(data$DPRK.condition %% 2 == 1, 1, 0)
# correct answer for the manipulation check on the DPRK question
data$DPRK.check.answer <- ifelse(data$DPRK.condition == 1 | data$DPRK.condition == 2, 1, 
                                 ifelse(data$DPRK.condition == 5 | data$DPRK.condition == 6, 2, 
                                        ifelse(data$DPRK.condition == 3 | data$DPRK.condition == 4, 3, 4)))

# female condition in the public assistance question
data$welfare.female <- ifelse(data$welfare.condition > 4, 1, 0)
# liberal condition in the public assistance question
data$welfare.liberal <- ifelse(data$welfare.condition == 3 | data$welfare.condition == 4 | data$welfare.condition == 7 | data$welfare.condition == 8, 1, 0)
# aggressuve condition in the public assistance question
data$welfare.aggressive <- ifelse(data$welfare.condition %% 2 == 1, 1, 0)
# correct answer for the manipulation check on the public assistance question
data$welfare.check.answer <- ifelse(data$welfare.condition == 1 | data$welfare.condition == 2, 1, 
                                    ifelse(data$welfare.condition == 5 | data$welfare.condition == 6, 2, 
                                           ifelse(data$welfare.condition == 3 | data$welfare.condition == 4, 3, 4)))

## manipulation checks
# correct answer for the manipulation check on the DPRK question
data$DPRK.correct <- ifelse((data$DPRK.female == 0 & data$DPRK.liberal == 0 & data$DPRK.check == 1) | 
                              (data$DPRK.female == 1 & data$DPRK.liberal == 0 & data$DPRK.check == 2) | 
                              (data$DPRK.female == 0 & data$DPRK.liberal == 1 & data$DPRK.check == 3) | 
                              (data$DPRK.female == 1 & data$DPRK.liberal == 1 & data$DPRK.check == 4), 1, 0)
# correct answer for the manipulation check on the public assistance question
data$welfare.correct <- ifelse((data$welfare.female == 0 & data$welfare.liberal == 0 & data$welfare.check == 1) | 
                                 (data$welfare.female == 1 & data$welfare.liberal == 0 & data$welfare.check == 2) | 
                                 (data$welfare.female == 0 & data$welfare.liberal == 1 & data$welfare.check == 3) | 
                                 (data$welfare.female == 1 & data$welfare.liberal == 1 & data$welfare.check == 4), 1, 0)

## respondents' ideology (Online Appendix, Table)
issue.attitudes <- data[, c("nuclear.power", "defense", "constitution", "public.security", 
                            "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax")]
attentive.issue.attitudes <- issue.attitudes[data$DPRK.correct == 1 | data$welfare.correct == 1 , ]
nrow(attentive.issue.attitudes)

# applying the ordinal IRT model to responses to the eight issue items
GRM.result <- grm(attentive.issue.attitudes)
# estimated discrimination parameters (Table A.9)
round(sapply(GRM.result$coefficients, function(x) x[5]), 2)
# respondents' ideological traits
GRM.score.result <- factor.scores(GRM.result, resp.patterns = issue.attitudes)
data$R.ideology <- GRM.score.result$score.dat$z1 * -1

## pooled data
pooled.data <- rbind(data.frame(outcome = data$DPRK.outcome, female = data$DPRK.female, 
                                liberal = data$DPRK.liberal, aggressive = data$DPRK.aggressive, 
                                order = data$order, welfare = 0, 
                                id = data$id, correct = data$DPRK.correct, R.gender = data$R.gender, 
                                age = data$age, middle.edu = data$middle.edu, high.edu = data$high.edu, 
                                R.ideology = data$R.ideology), 
                     data.frame(outcome = data$welfare.outcome, female = data$welfare.female, 
                                liberal = data$welfare.liberal, aggressive = data$welfare.aggressive, 
                                order = data$order, welfare = 1, 
                                id = data$id, correct = data$welfare.correct, R.gender = data$R.gender, 
                                age = data$age, middle.edu = data$middle.edu, high.edu = data$high.edu, 
                                R.ideology = data$R.ideology))

#### covariate balance checks using equivalene testing (Online Appendix C.1) ####
## function for equivalene testing
equivarence.test <- function(data.tested,  # data frame containing variables
                             treatment,  # treatment variable's name (1: treatment, 0: control)
                             covariates,  # covariates' names
                             epsilon  # equivalence limits
) {
  treatment.vec <- data.tested[, treatment]  # vector of the treatment variable
  treatment.group <- data.tested[treatment.vec == 1, ]  # treatment group's dara frame
  control.group <- data.tested[treatment.vec == 0, ]  # control group's dara frame
  results <- matrix(NA, length(covariates), 3)
  p.values <- rep(NA, length(covariates))
  for (i in 1:length(covariates)) {
    if (covariates[i] == treatment) next
    x.t <- treatment.group[, covariates[i]]  # i-th covariate in the treatment group
    x.c <- control.group[, covariates[i]]  # i-th covariate in the control group
    n.t <- sum(! is.na(x.t))  # omitting missing values
    n.c <- sum(! is.na(x.c))
    mean.t <- mean(x.t, na.rm = TRUE)  # mean in the treatment group
    mean.c <- mean(x.c, na.rm = TRUE)  # variance in the control group
    var.t <- var(x.t, na.rm = TRUE)  # mean in the treatment group
    var.c <- var(x.c, na.rm = TRUE)  # variance in the control group
    pooled.sd <- sqrt(((n.t - 1) * var.t + (n.c - 1) * var.c) / (n.t + n.c - 2))  # pooled standard deviation
    t.stat <- sqrt(n.t * n.c * (n.t + n.c - 2) / (n.t + n.c)) * (mean.t - mean.c) / 
      sqrt(sum((x.t - mean.t) ^ 2, na.rm = TRUE) + sum((x.c - mean.c) ^ 2, na.rm = TRUE))  # test statistic
    ncp <- (n.t * n.c * epsilon ^ 2) / (n.t + n.c)  # noncentrality parameter
    results[i, 1] <- mean.t - mean.c  # raw difference in means
    results[i, 2] <- (mean.t - mean.c) / pooled.sd  # standardized difference in means
    p.values[i] <- pf(abs(t.stat) ^ 2, 1, n.t + n.c - 2, ncp)  # p-value in equivalence testing
  }
  results[, 3] <- p.adjust(p.values, method = "BH")  # p-value adjustment by false discovery rate
  rownames(results) <- covariates
  colnames(results) <- c("raw.dif", "std.dif", "p.value")
  results
}

## results of covariate balance checks (Table A.3)
# DPRK question
round(cbind(equivarence.test(data, "DPRK.female", 
                             c("DPRK.female", "DPRK.liberal", "DPRK.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36), 
            equivarence.test(data, "DPRK.liberal", 
                             c("DPRK.female", "DPRK.liberal", "DPRK.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36), 
            equivarence.test(data, "DPRK.aggressive", 
                             c("DPRK.female", "DPRK.liberal", "DPRK.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36)), 3)
# public assistance question
round(cbind(equivarence.test(data, "welfare.female", 
                             c("welfare.female", "welfare.liberal", "welfare.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36), 
            equivarence.test(data, "welfare.liberal", 
                             c("welfare.female", "welfare.liberal", "welfare.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36), 
            equivarence.test(data, "welfare.aggressive", 
                             c("welfare.female", "welfare.liberal", "welfare.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36)), 3)

#### testing Hypothesis 1 by simple differences in means (Online Appendix C.2) ####
## bar plots (Figure A.1)
# DPRK question
DPRK.male.N.entire <- sum(data$DPRK.female == 0)  # number of observations
DPRK.female.N.entire <- sum(data$DPRK.female == 1)
DPRK.male.mean.entire <- mean(data$DPRK.outcome[data$DPRK.female == 0])  # difference in means
DPRK.female.mean.entire <- mean(data$DPRK.outcome[data$DPRK.female == 1])
DPRK.male.se.entire <- sd(data$DPRK.outcome[data$DPRK.female == 0]) / sqrt(DPRK.male.N.entire)  # standard error
DPRK.female.se.entire <- sd(data$DPRK.outcome[data$DPRK.female == 1]) / sqrt(DPRK.female.N.entire)
# welfare question
welfare.male.N.entire <- sum(data$welfare.female == 0)  # number of observations
welfare.female.N.entire <- sum(data$welfare.female == 1)
welfare.male.mean.entire <- mean(data$welfare.outcome[data$welfare.female == 0])  # difference in means
welfare.female.mean.entire <- mean(data$welfare.outcome[data$welfare.female == 1])
welfare.male.se.entire <- sd(data$welfare.outcome[data$welfare.female == 0]) / sqrt(welfare.male.N.entire)  # standard error
welfare.female.se.entire <- sd(data$welfare.outcome[data$welfare.female == 1]) / sqrt(welfare.female.N.entire)

# drawing the figure
cairo_pdf("Figure_A1.pdf", width = 5.9, height = 3.5, pointsize = 8)
layout(matrix(1:2, 1, 2))
par(mar = c(5, 5, 2.5, 1))
plot(NULL, NULL, type = "n", bty = "l", xlim = c(0.5, 2.5), ylim = c(1, 5), main = "DPRK", 
     xlab = "Experimental condition", ylab = "Disapproval of policy statements", xaxt = "n", yaxt = "n")
abline(h = 1:5, lty = 3, col = "gray")
polygon(c(0.7, 1.3, 1.3, 0.7), c(1, 1, DPRK.female.mean.entire, DPRK.female.mean.entire), col = gray(0.8))
polygon(c(1.7, 2.3, 2.3, 1.7), c(1, 1, DPRK.male.mean.entire, DPRK.male.mean.entire), col = gray(0.5))
arrows(1, DPRK.female.mean.entire + welfare.female.se.entire * qt(0.975, df = DPRK.female.N.entire - 1), 
       1, DPRK.female.mean.entire + welfare.female.se.entire * qt(0.025, df = DPRK.female.N.entire - 1), 
       length = 0.02, angle = 90, code = 3)
arrows(2, DPRK.male.mean.entire + welfare.male.se.entire * qt(0.975, df = DPRK.male.N.entire - 1), 
       2, DPRK.male.mean.entire + welfare.male.se.entire * qt(0.025, df = DPRK.male.N.entire - 1), 
       length = 0.02, angle = 90, code = 3)
axis(1, at = c(1, 2), labels = c("Female", "Male"))
axis(2, at = 1:5)
plot(NULL, NULL, type = "n", bty = "l", xlim = c(0.5, 2.5), ylim = c(1, 5), main = "public assistance", 
     xlab = "Experimental condition", ylab = "Disapproval of policy statements", xaxt = "n", yaxt = "n")
abline(h = 1:5, lty = 3, col = "gray")
polygon(c(0.7, 1.3, 1.3, 0.7), c(1, 1, welfare.female.mean.entire, welfare.female.mean.entire), col = gray(0.8))
polygon(c(1.7, 2.3, 2.3, 1.7), c(1, 1, welfare.male.mean.entire, welfare.male.mean.entire), col = gray(0.5))
arrows(1, welfare.female.mean.entire + welfare.female.se.entire * qt(0.975, df = welfare.female.N.entire - 1), 
       1, welfare.female.mean.entire + welfare.female.se.entire * qt(0.025, df = welfare.female.N.entire - 1), 
       length = 0.02, angle = 90, code = 3)
arrows(2, welfare.male.mean.entire + welfare.male.se.entire * qt(0.975, df = welfare.male.N.entire - 1), 
       2, welfare.male.mean.entire + welfare.male.se.entire * qt(0.025, df = welfare.male.N.entire - 1), 
       length = 0.02, angle = 90, code = 3)
axis(1, at = c(1, 2), labels = c("Female", "Male"))
axis(2, at = 1:5)
dev.off()

## Welch's t-tests
# DPRK question
round(DPRK.male.mean.entire - DPRK.female.mean.entire, 2)  # difference in means
DPRK.t.test.entire <- t.test(DPRK.outcome ~ DPRK.female, data)
round(DPRK.t.test.entire$p.value, 3)  # p-value

# welfare question
round(welfare.male.mean.entire - welfare.female.mean.entire, 2)  # difference in means
welfare.t.test.entire <- t.test(welfare.outcome ~ welfare.female, data)
round(welfare.t.test.entire$p.value, 3)  # p-value

#### testing Hypothesis 1 by linear models ####
## DPRK question
# linear regression
DPRK.H1.lm.entire <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + DPRK.aggressive + order, data)
# testing coefficients with the robust standard errors
DPRK.H1.result.entire <- coeftest(DPRK.H1.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.H1.result.entire[, 1:2], 2), round(DPRK.H1.result.entire[, 4], 3))
# number of observations
length(DPRK.H1.lm.entire$residuals)

## public assistance question
# linear regression
welfare.H1.lm.entire <- lm(welfare.outcome ~ welfare.female + welfare.liberal + welfare.aggressive + order, data)
# testing coefficients with the robust standard errors
welfare.H1.result.entire <- coeftest(welfare.H1.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.H1.result.entire[, 1:2], 2), round(welfare.H1.result.entire[, 4], 3))
# number of observations
length(welfare.H1.lm.entire$residuals)

## pooled analysis
# linear regression
pooled.H1.lm.entire <- lm(outcome ~ female + liberal + aggressive + order + welfare, pooled.data)
# testing coefficients with the robust standard errors
pooled.H1.result.entire <- coef_test(pooled.H1.lm.entire, 
                                     vcov = vcovCR(pooled.H1.lm.entire, 
                                                   cluster = pooled.data[, "id"], 
                                                   type = "CR2"))
# summary
cbind(round(pooled.H1.result.entire[, 1:2], 2), p = round(pooled.H1.result.entire$p_Satt, 3))
# number of observations
length(pooled.H1.lm.entire$residuals)

#### testing Hypothesis 3 ####
## DPRK question
# linear regression
DPRK.H3.lm.entire <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + DPRK.aggressive + order + DPRK.female : DPRK.liberal, data)
# testing coefficients with the robust standard errors
DPRK.H3.result.entire <- coeftest(DPRK.H3.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.H3.result.entire[, 1:2], 2), round(DPRK.H3.result.entire[, 4], 3))
# number of observations
length(DPRK.H3.lm.entire$residuals)

## welfare question
# linear regression
welfare.H3.lm.entire <- lm(welfare.outcome ~ welfare.female + welfare.liberal + welfare.aggressive + order + welfare.female : welfare.liberal, data)
# testing coefficients with the robust standard errors
welfare.H3.result.entire <- coeftest(welfare.H3.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.H3.result.entire[, 1:2], 2), round(welfare.H3.result.entire[, 4], 3))
# number of observations
length(welfare.H3.lm.entire$residuals)

## pooled analysis
# linear regression
pooled.H3.lm.entire <- lm(outcome ~ female + liberal + aggressive + order + female : liberal + welfare, pooled.data)
# testing coefficients with the robust standard errors
pooled.H3.result.entire <- coef_test(pooled.H3.lm.entire, 
                                     vcov = vcovCR(pooled.H3.lm.entire, 
                                                   cluster = pooled.data[, "id"], 
                                                   type = "CR2"))
# summary
cbind(round(pooled.H3.result.entire[, 1:2], 2), p = round(pooled.H3.result.entire$p_Satt, 3))
# number of observations
length(pooled.H3.lm.entire$residuals)

#### summarizing the results of testing Hypotheses 1 and 3 (Figure 1) ####
## confidence intervals for the tests of Hypothesis 1
DPRK.H1.CI.entire <- coefci(DPRK.H1.lm.entire, vcov. = vcovHC, type = "HC2")
welfare.H1.CI.entire <- coefci(welfare.H1.lm.entire, vcov. = vcovHC, type = "HC2")
pooled.H1.CI.entire <- conf_int(pooled.H1.lm.entire, 
                                vcov = vcovCR(pooled.H1.lm.entire, 
                                              cluster = pooled.data[, "id"], 
                                              type = "CR2"))

## confidence intervals for the tests of Hypothesis 3
# effect of candidate gender for the conservative condition
DPRK.H3.conservative.CI.entire <- coefci(DPRK.H3.lm.entire, vcov. = vcovHC, type = "HC2")
welfare.H3.conservative.CI.entire <- coefci(welfare.H3.lm.entire, vcov. = vcovHC, type = "HC2")
pooled.H3.conservative.CI.entire <- conf_int(pooled.H3.lm.entire, 
                                             vcov = vcovCR(pooled.H3.lm.entire, 
                                                           cluster = pooled.data[, "id"], 
                                                           type = "CR2"))

# effect of candidate gender for the liberal condition
DPRK.H3.lm.entire.reversed <- lm(DPRK.outcome ~ DPRK.female + I(1 - DPRK.liberal) + DPRK.aggressive + 
                                    order + DPRK.female : I(1 - DPRK.liberal), data)
DPRK.H3.liberal.CI.entire <- coefci(DPRK.H3.lm.entire.reversed, vcov. = vcovHC, type = "HC2")
welfare.H3.lm.entire.reversed <- lm(welfare.outcome ~ welfare.female + I(1 - welfare.liberal) + welfare.aggressive + 
                                       order + welfare.female : I(1 - welfare.liberal), data)
welfare.H3.liberal.CI.entire <- coefci(welfare.H3.lm.entire.reversed, vcov. = vcovHC, type = "HC2")
pooled.H3.lm.entire.reversed <- lm(outcome ~ female + I(1 - liberal) + aggressive + order + 
                                      female : I(1 - liberal) + welfare, pooled.data)
pooled.H3.liberal.CI.entire <- conf_int(pooled.H3.lm.entire.reversed, 
                                         vcov = vcovCR(pooled.H3.lm.entire.reversed, 
                                                       cluster = pooled.data[, "id"], 
                                                       type = "CR2"))

png("Figure_1.png", width = 4.7, height = 3, units = "in", pointsize = 8, res = 1200)
par(mar = c(5, 1, 1, 1))
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.8, 0.7), ylim = c(0, 14), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = 0, lty = 3, col = "gray")
points(DPRK.H1.lm.entire$coefficients[2], 13, pch = 19)
segments(DPRK.H1.CI.entire[2, 1], 13, DPRK.H1.CI.entire[2, 2], 13)
points(DPRK.H3.lm.entire.reversed$coefficients[2], 12, pch = 19)
segments(DPRK.H3.liberal.CI.entire[2, 1], 12, DPRK.H3.liberal.CI.entire[2, 2], 12)
points(DPRK.H3.lm.entire$coefficients[2], 11, pch = 19)
segments(DPRK.H3.conservative.CI.entire[2, 1], 11, DPRK.H3.conservative.CI.entire[2, 2], 11)
segments(DPRK.H3.liberal.CI.entire[2, 2] + 0.025, 12, DPRK.H3.liberal.CI.entire[2, 2] + 0.1, 12, lwd = 0.5)
segments(DPRK.H3.conservative.CI.entire[2, 2] + 0.025, 11, DPRK.H3.liberal.CI.entire[2, 2] + 0.1, 11, lwd = 0.5)
segments(DPRK.H3.liberal.CI.entire[2, 2] + 0.1, 12, DPRK.H3.liberal.CI.entire[2, 2] + 0.1, 11, lwd = 0.5)
text(DPRK.H3.liberal.CI.entire[2, 2] + 0.1, 11.5, 
     paste0(sprintf("%4.2f", DPRK.H3.lm.entire$coefficients[6]), " [\u2212", 
            sprintf("%4.2f", abs(DPRK.H3.conservative.CI.entire[6, 1])), ", ", 
            sprintf("%4.2f", DPRK.H3.conservative.CI.entire[6, 2]), "]"), cex = 0.8, pos = 4)
points(welfare.H1.lm.entire$coefficients[2], 8, pch = 19)
segments(welfare.H1.CI.entire[2, 1], 8, welfare.H1.CI.entire[2, 2], 8)
points(welfare.H3.lm.entire.reversed$coefficients[2], 7, pch = 19)
segments(welfare.H3.liberal.CI.entire[2, 1], 7, welfare.H3.liberal.CI.entire[2, 2], 7)
points(welfare.H3.lm.entire$coefficients[2], 6, pch = 19)
segments(welfare.H3.conservative.CI.entire[2, 1], 6, welfare.H3.conservative.CI.entire[2, 2], 6)
segments(welfare.H3.liberal.CI.entire[2, 2] + 0.025, 7, welfare.H3.conservative.CI.entire[2, 2] + 0.1, 7, lwd = 0.5)
segments(welfare.H3.conservative.CI.entire[2, 2] + 0.025, 6, welfare.H3.conservative.CI.entire[2, 2] + 0.1, 6, lwd = 0.5)
segments(welfare.H3.conservative.CI.entire[2, 2] + 0.1, 7, welfare.H3.conservative.CI.entire[2, 2] + 0.1, 6, lwd = 0.5)
text(welfare.H3.conservative.CI.entire[2, 2] + 0.1, 6.5, 
     paste0("\u2212", sprintf("%4.2f", abs(welfare.H3.lm.entire$coefficients[6])), " [\u2212", 
            sprintf("%4.2f", abs(welfare.H3.conservative.CI.entire[6, 1])), ", ", 
            sprintf("%4.2f", welfare.H3.conservative.CI.entire[6, 2]), "]"), cex = 0.8, pos = 4)
points(pooled.H1.lm.entire$coefficients[2], 3, pch = 19)
segments(pooled.H1.CI.entire[2, 4], 3, pooled.H1.CI.entire[2, 5], 3)
points(pooled.H3.lm.entire.reversed$coefficients[2], 2, pch = 19)
segments(pooled.H3.liberal.CI.entire[2, 4], 2, pooled.H3.liberal.CI.entire[2, 5], 2)
points(pooled.H3.lm.entire$coefficients[2], 1, pch = 19)
segments(pooled.H3.conservative.CI.entire[2, 4], 1, pooled.H3.conservative.CI.entire[2, 5], 1)
segments(pooled.H3.liberal.CI.entire[2, 5] + 0.025, 2, pooled.H3.conservative.CI.entire[2, 5] + 0.1, 2, lwd = 0.5)
segments(pooled.H3.conservative.CI.entire[2, 5] + 0.025, 1, pooled.H3.conservative.CI.entire[2, 5] + 0.1, 1, lwd = 0.5)
segments(pooled.H3.conservative.CI.entire[2, 5] + 0.1, 2, pooled.H3.conservative.CI.entire[2, 5] + 0.1, 1, lwd = 0.5)
text(pooled.H3.conservative.CI.entire[2, 5] + 0.1, 1.5, 
     paste0(sprintf("%4.2f", abs(pooled.H3.lm.entire$coefficients[7])), " [\u2212", 
            sprintf("%4.2f", abs(pooled.H3.conservative.CI.entire[7, 4])), ", ", 
            sprintf("%4.2f", pooled.H3.conservative.CI.entire[7, 5]), "]"), cex = 0.8, pos = 4)
text(rep(-0.8, 3), c(14, 9, 4), c("DPRK", "Public assistance", "Pooled"), pos = 4)
text(rep(-0.75, 3), c(13, 8, 3), rep("Female"), pos = 4)
text(rep(-0.75, 3), c(12, 7, 2), rep("Female | liberal"), pos = 4)
text(rep(-0.75, 3), c(11, 6, 1), rep("Female | conservative"), pos = 4)
axis(1, at = seq(-0.2, 0.6, 0.2))
mtext("Effects", side = 1, at = 0.2, line = 3)
dev.off()

#### testing heterogenous treatment effects by respondent gender (Online Appendix C.3) ####
## DPRK question
# linear regression
DPRK.hetero.gender.lm.entire <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + 
                                     DPRK.aggressive + order + R.gender + R.gender : DPRK.female, data)
# testing coefficients with the robust standard errors
DPRK.hetero.gender.result.entire <- coeftest(DPRK.hetero.gender.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.hetero.gender.result.entire[, 1:2], 2), round(DPRK.hetero.gender.result.entire[, 4], 3))
# number of observations
length(DPRK.hetero.gender.lm.entire$residuals)

## public assistance question
# linear regression
welfare.hetero.gender.lm.entire <- lm(welfare.outcome ~ welfare.female + welfare.liberal + 
                                        welfare.aggressive + order + R.gender + R.gender : welfare.female, data)
# testing coefficients with the robust standard errors
welfare.hetero.gender.result.entire <- coeftest(welfare.hetero.gender.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.hetero.gender.result.entire[, 1:2], 2), round(welfare.hetero.gender.result.entire[, 4], 3))
# number of observations
length(welfare.hetero.gender.lm.entire$residuals)

## pooled analysis
# linear regression
pooled.hetero.gender.lm.entire <- lm(outcome ~ female + liberal + aggressive + order + welfare + 
                                       R.gender + R.gender : female, pooled.data)
# testing coefficients with the robust standard errors
pooled.hetero.gender.result.entire <- coef_test(pooled.hetero.gender.lm.entire, 
                                                vcov = vcovCR(pooled.hetero.gender.lm.entire, 
                                                              cluster = pooled.data[, "id"], 
                                                              type = "CR2"))
# summary
cbind(round(pooled.hetero.gender.result.entire[, 1:2], 2), p = round(pooled.hetero.gender.result.entire$p_Satt, 3))
# number of observations
length(pooled.hetero.gender.lm.entire$residuals)

#### testing heterogenous treatment effects by respondent ideology (Online Appendix C.4) ####
# standardize respondents' ideology
data$R.ideology.std.entire <- (data$R.ideology - mean(data$R.ideology)) / sd(data$R.ideology)

## DPRK question
# linear regression
DPRK.hetero.ideology.lm.entire <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + 
                                       DPRK.aggressive + order + R.ideology.std.entire + R.ideology.std.entire : DPRK.female, data)
# testing coefficients with the robust standard errors
DPRK.hetero.ideology.result.entire <- coeftest(DPRK.hetero.ideology.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.hetero.ideology.result.entire[, 1:2], 2), round(DPRK.hetero.ideology.result.entire[, 4], 3))
# number of observations
length(DPRK.hetero.ideology.lm.entire$residuals)

## public assistance question
# linear regression
welfare.hetero.ideology.lm.entire <- lm(welfare.outcome ~ welfare.female + welfare.liberal + 
                                          welfare.aggressive + order + R.ideology.std.entire + R.ideology.std.entire : welfare.female, data)
# testing coefficients with the robust standard errors
welfare.hetero.ideology.result.entire <- coeftest(welfare.hetero.ideology.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.hetero.ideology.result.entire[, 1:2], 2), round(welfare.hetero.ideology.result.entire[, 4], 3))
# number of observations
length(welfare.hetero.ideology.lm.entire$residuals)

## pooled analysis
# standardize respondents' ideology
pooled.data$R.ideology.std.entire <- (pooled.data$R.ideology - mean(pooled.data$R.ideology)) / sd(pooled.data$R.ideology)
# linear regression
pooled.hetero.ideology.lm.entire <- lm(outcome ~ female + liberal + aggressive + order + welfare + 
                                         R.ideology.std.entire + R.ideology.std.entire : female, pooled.data)
# testing coefficients with the robust standard errors
pooled.hetero.ideology.result.entire <- coef_test(pooled.hetero.ideology.lm.entire, 
                                                  vcov = vcovCR(pooled.hetero.ideology.lm.entire, 
                                                                cluster = pooled.data[, "id"], 
                                                                type = "CR2"))
# summary
cbind(round(pooled.hetero.ideology.result.entire[, 1:2], 2), p = round(pooled.hetero.ideology.result.entire$p_Satt, 3))
# number of observations
length(pooled.hetero.ideology.lm.entire$residuals)

#### testing Hypothesis 4 (Online Appendix C.5) ####
## DPRK question
# linear regression
DPRK.H4.lm.entire <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + DPRK.aggressive + order + DPRK.female : DPRK.aggressive, data)
# testing coefficients with the robust standard errors
DPRK.H4.result.entire <- coeftest(DPRK.H4.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.H4.result.entire[, 1:2], 2), round(DPRK.H4.result.entire[, 4], 3))
# number of observations
length(DPRK.H4.lm.entire$residuals)

## welfare question
# linear regression
welfare.H4.lm.entire <- lm(welfare.outcome ~ welfare.female + welfare.liberal + welfare.aggressive + order + welfare.female : welfare.aggressive, data)
# testing coefficients with the robust standard errors
welfare.H4.result.entire <- coeftest(welfare.H4.lm.entire, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.H4.result.entire[, 1:2], 2), round(welfare.H4.result.entire[, 4], 3))
# number of observations
length(welfare.H4.lm.entire$residuals)

## pooled analysis
# linear regression
pooled.H4.lm.entire <- lm(outcome ~ female + liberal + aggressive + order + female : aggressive + welfare, pooled.data)
# testing coefficients with the robust standard errors
pooled.H4.result.entire <- coef_test(pooled.H4.lm.entire, 
                                      vcov = vcovCR(pooled.H4.lm.entire, 
                                                    cluster = pooled.data[, "id"], 
                                                    type = "CR2"))
# summary
cbind(round(pooled.H4.result.entire[, 1:2], 2), p = round(pooled.H4.result.entire$p_Satt, 3))
# number of observations
length(pooled.H4.lm.entire$residuals)

#### details of the factual manipulation checks (Online Appendix D) ####
## relationship between the group assignment and responses for the manipulation checks (Table A.8)
# DPRK question
table(data$DPRK.check.answer, data$DPRK.check)
round(t(apply(table(data$DPRK.check.answer, data$DPRK.check), 1, prop.table)), 3) * 100
table(data$DPRK.check)
round(prop.table(table(data$DPRK.check)), 3) * 100
table(data$DPRK.check.answer)
sum(table(data$DPRK.check.answer))
# welfare question
table(data$welfare.check.answer, data$welfare.check)
round(t(apply(table(data$welfare.check.answer, data$welfare.check), 1, prop.table)), 3) * 100
table(data$welfare.check)
round(prop.table(table(data$welfare.check)), 3) * 100
table(data$welfare.check.answer)
sum(table(data$welfare.check.answer))

## results of covariate balance checks without satisficers (Table A.9)
round(cbind(equivarence.test(data[data$DPRK.correct == 1, ], "DPRK.female", 
                             c("DPRK.female", "DPRK.liberal", "DPRK.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36), 
            equivarence.test(data[data$DPRK.correct == 1, ], "DPRK.liberal", 
                             c("DPRK.female", "DPRK.liberal", "DPRK.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36), 
            equivarence.test(data[data$DPRK.correct == 1, ], "DPRK.aggressive", 
                             c("DPRK.female", "DPRK.liberal", "DPRK.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36)), 3)

round(cbind(equivarence.test(data[data$welfare.correct == 1, ], "welfare.female", 
                             c("welfare.female", "welfare.liberal", "welfare.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36), 
            equivarence.test(data[data$welfare.correct == 1, ], "welfare.liberal", 
                             c("welfare.female", "welfare.liberal", "welfare.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36), 
            equivarence.test(data[data$welfare.correct == 1, ], "welfare.aggressive", 
                             c("welfare.female", "welfare.liberal", "welfare.aggressive", "order", 
                               "R.gender", "age", "middle.edu", "high.edu", "knowledge", 
                               "nuclear.power", "defense", "constitution", "public.security", 
                               "separate.surnames", "preemptive.attack", "same.sex.marriage", "income.tax"), 
                             0.36)), 3)

## distribution of outcomes of attentive and inattentive respondents (Fugire A.2)
# setting colors
Paired <- colorschemes$Categorical.12

# drawing the figure
cairo_pdf("Figure_A2.pdf", width = 6, height = 3, pointsize = 7)
layout(matrix(1:2, 1, 2, byrow = TRUE))
par(mar = c(4, 4, 3, 1), lwd = 0.5)
hist(data$DPRK.outcome[data$DPRK.correct == 0], breaks = 1:6, freq = FALSE, right = FALSE, 
     col = paste0(Paired[11], "80"), border = Paired[12], density = 50, ylim = c(0, 0.5), axes = FALSE, 
     main = "DPRK question", xlab = "", ylab = "Proportion")
hist(data$DPRK.outcome[data$DPRK.correct == 1], breaks = 1:6, freq = FALSE, right = FALSE, 
     col = paste0(Paired[9], "80"), border = Paired[10], add = TRUE)
legend("topright", legend = c("Attentive", "Inattentive"), 
       fill = c(paste0(Paired[9], "80"), paste0(Paired[11], "80")), 
       border = c(Paired[10], Paired[12]), bty = "n", density = c(NA, 50))
segments(1, 0, 5, 0)
axis(2, lwd = 0.5)
mtext(c("Approve\n\n", "Somewhat\napprove\n", "Neither\napprove nor\ndisapprove", 
        "Somewhat\ndisapprove\n", "Disapprove\n\n"), side = 1, at = 1:5 + 0.5, line = 1, cex = 0.65)
mtext("Outcome", side = 1, line = 2.5)
hist(data$welfare.outcome[data$welfare.correct == 0], breaks = 1:6, freq = FALSE, right = FALSE, 
     col = paste0(Paired[11], "80"), border = Paired[12], density = 50, ylim = c(0, 0.5), axes = FALSE, 
     main = "public assistance question", xlab = "", ylab = "Proportion")
hist(data$welfare.outcome[data$welfare.correct == 1], breaks = 1:6, freq = FALSE, right = FALSE, 
     col = paste0(Paired[9], "80"), border = Paired[10], add = TRUE)
legend("topright", legend = c("Attentive", "Inattentive"), 
       fill = c(paste0(Paired[9], "80"), paste0(Paired[11], "80")), 
       border = c(Paired[10], Paired[12]), bty = "n", density = c(NA, 50))
segments(1, 0, 5, 0)
axis(2, lwd = 0.5)
mtext(c("Approve\n\n", "Somewhat\napprove\n", "Neither\napprove nor\ndisapprove", 
        "Somewhat\ndisapprove\n", "Disapprove\n\n"), side = 1, at = 1:5 + 0.5, line = 1, cex = 0.65)
mtext("Outcome", side = 1, line = 2.5)
dev.off()

#### Analyses using the data from attentive respondents (Online Appendix E) ####
### testing Hypothesis 1 by simple differences in means
## bar plots
# DPRK question
DPRK.male.N.partial <- sum(data$DPRK.female == 0 & data$DPRK.correct == 1)  # number of observations
DPRK.female.N.partial <- sum(data$DPRK.female == 1 & data$DPRK.correct == 1)
DPRK.male.mean.partial <- mean(data$DPRK.outcome[data$DPRK.female == 0 & data$DPRK.correct == 1])  # difference in means
DPRK.female.mean.partial <- mean(data$DPRK.outcome[data$DPRK.female == 1 & data$DPRK.correct == 1])
DPRK.male.se.partial <- sd(data$DPRK.outcome[data$DPRK.female == 0 & data$DPRK.correct == 1]) / sqrt(DPRK.male.N.partial)  # standard error
DPRK.female.se.partial <- sd(data$DPRK.outcome[data$DPRK.female == 1 & data$DPRK.correct == 1]) / sqrt(DPRK.female.N.partial)
# welfare question
welfare.male.N.partial <- sum(data$welfare.female == 0 & data$welfare.correct == 1)  # number of observations
welfare.female.N.partial <- sum(data$welfare.female == 1 & data$welfare.correct == 1)
welfare.male.mean.partial <- mean(data$welfare.outcome[data$welfare.female == 0 & data$welfare.correct == 1])  # difference in means
welfare.female.mean.partial <- mean(data$welfare.outcome[data$welfare.female == 1 & data$welfare.correct == 1])
welfare.male.se.partial <- sd(data$welfare.outcome[data$welfare.female == 0 & data$welfare.correct == 1]) / sqrt(welfare.male.N.partial)  # standard error
welfare.female.se.partial <- sd(data$welfare.outcome[data$welfare.female == 1 & data$welfare.correct == 1]) / sqrt(welfare.female.N.partial)

# drawing the figure
cairo_pdf("Figure_A3.pdf", width = 5.9, height = 3.5, pointsize = 8)
layout(matrix(1:2, 1, 2))
par(mar = c(5, 5, 2.5, 1))
plot(NULL, NULL, type = "n", bty = "l", xlim = c(0.5, 2.5), ylim = c(1, 5), main = "DPRK", 
     xlab = "Experimental condition", ylab = "Disapproval of policy statements", xaxt = "n", yaxt = "n")
abline(h = 1:5, lty = 3, col = "gray")
polygon(c(0.7, 1.3, 1.3, 0.7), c(1, 1, DPRK.female.mean.partial, DPRK.female.mean.partial), col = gray(0.8))
polygon(c(1.7, 2.3, 2.3, 1.7), c(1, 1, DPRK.male.mean.partial, DPRK.male.mean.partial), col = gray(0.5))
arrows(1, DPRK.female.mean.partial + welfare.female.se.partial * qt(0.975, df = DPRK.female.N.partial - 1), 
       1, DPRK.female.mean.partial + welfare.female.se.partial * qt(0.025, df = DPRK.female.N.partial - 1), 
       length = 0.02, angle = 90, code = 3)
arrows(2, DPRK.male.mean.partial + welfare.male.se.partial * qt(0.975, df = DPRK.male.N.partial - 1), 
       2, DPRK.male.mean.partial + welfare.male.se.partial * qt(0.025, df = DPRK.male.N.partial - 1), 
       length = 0.02, angle = 90, code = 3)
axis(1, at = c(1, 2), labels = c("Female", "Male"))
axis(2, at = 1:5)
plot(NULL, NULL, type = "n", bty = "l", xlim = c(0.5, 2.5), ylim = c(1, 5), main = "public assistance", 
     xlab = "Experimental condition", ylab = "Disapproval of policy statements", xaxt = "n", yaxt = "n")
abline(h = 1:5, lty = 3, col = "gray")
polygon(c(0.7, 1.3, 1.3, 0.7), c(1, 1, welfare.female.mean.partial, welfare.female.mean.partial), col = gray(0.8))
polygon(c(1.7, 2.3, 2.3, 1.7), c(1, 1, welfare.male.mean.partial, welfare.male.mean.partial), col = gray(0.5))
arrows(1, welfare.female.mean.partial + welfare.female.se.partial * qt(0.975, df = welfare.female.N.partial - 1), 
       1, welfare.female.mean.partial + welfare.female.se.partial * qt(0.025, df = welfare.female.N.partial - 1), 
       length = 0.02, angle = 90, code = 3)
arrows(2, welfare.male.mean.partial + welfare.male.se.partial * qt(0.975, df = welfare.male.N.partial - 1), 
       2, welfare.male.mean.partial + welfare.male.se.partial * qt(0.025, df = welfare.male.N.partial - 1), 
       length = 0.02, angle = 90, code = 3)
axis(1, at = c(1, 2), labels = c("Female", "Male"))
axis(2, at = 1:5)
dev.off()

## Welch's t-tests
# DPRK question
round(DPRK.male.mean.partial - DPRK.female.mean.partial, 2)  # difference in means
DPRK.t.test.partial <- t.test(DPRK.outcome ~ DPRK.female, data, subset = DPRK.correct == 1)
round(DPRK.t.test.partial$p.value, 3)  # p-value

# welfare question
round(welfare.male.mean.partial - welfare.female.mean.partial, 2)  # difference in means
welfare.t.test.partial <- t.test(welfare.outcome ~ welfare.female, data, subset = welfare.correct == 1)
round(welfare.t.test.partial$p.value, 3)  # p-value

## Cohen's d for the public assistance question
welfare.cohen.d.partial <- cohen.d(data$welfare.outcome[data$welfare.correct == 1], 
                                   factor(data$welfare.female[data$welfare.correct == 1]))
round(welfare.cohen.d.partial$estimate, 2)  # Cohen's d
round(welfare.cohen.d.partial$sd, 2)  # pooled standard deviation

### testing Hypothesis 1 by linear models
## DPRK question
# linear regression
DPRK.H1.lm.partial <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + DPRK.aggressive + order, 
                         data, subset = DPRK.correct == 1)
# testing coefficients with the robust standard errors
DPRK.H1.result.partial <- coeftest(DPRK.H1.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.H1.result.partial[, 1:2], 2), round(DPRK.H1.result.partial[, 4], 3))
# number of observations
length(DPRK.H1.lm.partial$residuals)

## public assistance question
# linear regression
welfare.H1.lm.partial <- lm(welfare.outcome ~ welfare.female + welfare.liberal + welfare.aggressive + order, 
                            data, subset = welfare.correct == 1)
# testing coefficients with the robust standard errors
welfare.H1.result.partial <- coeftest(welfare.H1.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.H1.result.partial[, 1:2], 2), round(welfare.H1.result.partial[, 4], 3))
# number of observations
length(welfare.H1.lm.partial$residuals)

## pooled analysis
# linear regression
pooled.H1.lm.partial <- lm(outcome ~ female + liberal + aggressive + order + welfare, 
                           pooled.data, subset = correct == 1)
# testing coefficients with the robust standard errors
pooled.H1.result.partial <- coef_test(pooled.H1.lm.partial, 
                                      vcov = vcovCR(pooled.H1.lm.partial, 
                                                    cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                    type = "CR2"))
# summary
cbind(round(pooled.H1.result.partial[, 1:2], 2), p = round(pooled.H1.result.partial$p_Satt, 3))
# number of observations
length(pooled.H1.lm.partial$residuals)

### testing Hypothesis 2
# linear regression
pooled.H2.lm.partial <- lm(outcome ~ female + liberal + aggressive + order + welfare + female : welfare, 
                           pooled.data, subset = correct == 1)
# testing coefficients with the robust standard errors
pooled.H2.result.partial <- coef_test(pooled.H2.lm.partial, 
                                      vcov = vcovCR(pooled.H2.lm.partial, 
                                                    cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                    type = "CR2"))
# summary
cbind(round(pooled.H2.result.partial[, 1:2], 2), p = round(pooled.H2.result.partial$p_Satt, 3))
# number of observations
length(pooled.H2.lm.partial$residuals)

# conditional effect
pooled.H2.lm.partial.reversed <- lm(outcome ~ female + liberal + aggressive + order + 
                                      I(1 - welfare) + female : I(1 - welfare), 
                                    pooled.data, subset = correct == 1)
pooled.H2.result.partial.reversed <- coef_test(pooled.H2.lm.partial.reversed, 
                                               vcov = vcovCR(pooled.H2.lm.partial.reversed, 
                                                             cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                             type = "CR2"))
cbind(round(pooled.H2.result.partial.reversed[, 1:2], 2), p = round(pooled.H2.result.partial.reversed$p_Satt, 3))

### testing heterogenous treatment effects by respondent gender
## DPRK question
# linear regression
DPRK.hetero.gender.lm.partial <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + 
                                      DPRK.aggressive + order + R.gender + R.gender : DPRK.female, 
                                    data, subset = DPRK.correct == 1)
# testing coefficients with the robust standard errors
DPRK.hetero.gender.result.partial <- coeftest(DPRK.hetero.gender.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.hetero.gender.result.partial[, 1:2], 2), round(DPRK.hetero.gender.result.partial[, 4], 3))
# number of observations
length(DPRK.hetero.gender.lm.partial$residuals)

## public assistance question
# linear regression
welfare.hetero.gender.lm.partial <- lm(welfare.outcome ~ welfare.female + welfare.liberal + 
                                         welfare.aggressive + order + R.gender + R.gender : welfare.female, 
                                       data, subset = welfare.correct == 1)
# testing coefficients with the robust standard errors
welfare.hetero.gender.result.partial <- coeftest(welfare.hetero.gender.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.hetero.gender.result.partial[, 1:2], 2), round(welfare.hetero.gender.result.partial[, 4], 3))
# number of observations
length(welfare.hetero.gender.lm.partial$residuals)

## pooled analysis
# linear regression
pooled.hetero.gender.lm.partial <- lm(outcome ~ female + liberal + aggressive + order + welfare + 
                                        R.gender + R.gender : female, 
                                      pooled.data, subset = correct == 1)
# testing coefficients with the robust standard errors
pooled.hetero.gender.result.partial <- coef_test(pooled.hetero.gender.lm.partial, 
                                                 vcov = vcovCR(pooled.hetero.gender.lm.partial, 
                                                               cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                               type = "CR2"))
# summary
cbind(round(pooled.hetero.gender.result.partial[, 1:2], 2), p = round(pooled.hetero.gender.result.partial$p_Satt, 3))
# number of observations
length(pooled.hetero.gender.lm.partial$residuals)

### testing heterogenous treatment effects by respondent ideology
## DPRK question
# standardize respondents' ideology
data$R.ideology.std.partial.DPRK <- (data$R.ideology - mean(data$R.ideology[data$DPRK.correct == 1])) / 
  sd(data$R.ideology[data$DPRK.correct == 1])
# linear regression
DPRK.hetero.ideology.lm.partial <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + 
                                        DPRK.aggressive + order + R.ideology.std.partial.DPRK + R.ideology.std.partial.DPRK : DPRK.female, 
                                      data, subset = DPRK.correct == 1)
# testing coefficients with the robust standard errors
DPRK.hetero.ideology.result.partial <- coeftest(DPRK.hetero.ideology.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.hetero.ideology.result.partial[, 1:2], 2), round(DPRK.hetero.ideology.result.partial[, 4], 3))
# number of observations
length(DPRK.hetero.ideology.lm.partial$residuals)

## public assistance question
# standardize respondents' ideology
data$R.ideology.std.partial.welfare <- (data$R.ideology - mean(data$R.ideology[data$welfare.correct == 1])) / 
  sd(data$R.ideology[data$welfare.correct == 1])
# linear regression
welfare.hetero.ideology.lm.partial <- lm(welfare.outcome ~ welfare.female + welfare.liberal + 
                                           welfare.aggressive + order + R.ideology.std.partial.welfare + 
                                           R.ideology.std.partial.welfare : welfare.female, 
                                         data, subset = welfare.correct == 1)
# testing coefficients with the robust standard errors
welfare.hetero.ideology.result.partial <- coeftest(welfare.hetero.ideology.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.hetero.ideology.result.partial[, 1:2], 2), round(welfare.hetero.ideology.result.partial[, 4], 3))
# number of observations
length(welfare.hetero.ideology.lm.partial$residuals)

## pooled analysis
# standardize respondents' ideology
pooled.data$R.ideology.std.partial <- (pooled.data$R.ideology - mean(pooled.data$R.ideology[pooled.data$correct == 1])) / 
  sd(pooled.data$R.ideology[pooled.data$correct == 1])
# linear regression
pooled.hetero.ideology.lm.partial <- lm(outcome ~ female + liberal + aggressive + order + welfare + 
                                          R.ideology.std.partial + R.ideology.std.partial : female, 
                                        pooled.data, subset = correct == 1)
# testing coefficients with the robust standard errors
pooled.hetero.ideology.result.partial <- coef_test(pooled.hetero.ideology.lm.partial, 
                                                   vcov = vcovCR(pooled.hetero.ideology.lm.partial, 
                                                                 cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                                 type = "CR2"))
# summary
cbind(round(pooled.hetero.ideology.result.partial[, 1:2], 2), p = round(pooled.hetero.ideology.result.partial$p_Satt, 3))
# number of observations
length(pooled.hetero.ideology.lm.partial$residuals)

### testing Hypothesis 3
## DPRK question
# linear regression
DPRK.H3.lm.partial <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + DPRK.aggressive + order + DPRK.female : DPRK.liberal, 
                         data, subset = DPRK.correct == 1)
# testing coefficients with the robust standard errors
DPRK.H3.result.partial <- coeftest(DPRK.H3.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.H3.result.partial[, 1:2], 2), round(DPRK.H3.result.partial[, 4], 3))
# number of observations
length(DPRK.H3.lm.partial$residuals)

## welfare question
# linear regression
welfare.H3.lm.partial <- lm(welfare.outcome ~ welfare.female + welfare.liberal + welfare.aggressive + order + welfare.female : welfare.liberal, 
                            data, subset = welfare.correct == 1)
# testing coefficients with the robust standard errors
welfare.H3.result.partial <- coeftest(welfare.H3.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.H3.result.partial[, 1:2], 2), round(welfare.H3.result.partial[, 4], 3))
# number of observations
length(welfare.H3.lm.partial$residuals)

## pooled analysis
# linear regression
pooled.H3.lm.partial <- lm(outcome ~ female + liberal + aggressive + order + female : liberal + welfare, 
                           pooled.data, subset = correct == 1)
# testing coefficients with the robust standard errors
pooled.H3.result.partial <- coef_test(pooled.H3.lm.partial, 
                                      vcov = vcovCR(pooled.H3.lm.partial, 
                                                    cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                    type = "CR2"))
# summary
cbind(round(pooled.H3.result.partial[, 1:2], 2), p = round(pooled.H3.result.partial$p_Satt, 3))
# number of observations
length(pooled.H3.lm.partial$residuals)

### summarizing the results of testing Hypotheses 1 and 3 (Figure A.4)
## confidence intervals for the tests of Hypothesis 1
DPRK.H1.CI.partial <- coefci(DPRK.H1.lm.partial, vcov. = vcovHC, type = "HC2")
welfare.H1.CI.partial <- coefci(welfare.H1.lm.partial, vcov. = vcovHC, type = "HC2")
pooled.H1.CI.partial <- conf_int(pooled.H1.lm.partial, 
                                 vcov = vcovCR(pooled.H1.lm.partial, 
                                               cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                               type = "CR2"))

## confidence intervals for the tests of Hypothesis 3
# effect of candidate gender for the conservative condition
DPRK.H3.conservative.CI.partial <- coefci(DPRK.H3.lm.partial, vcov. = vcovHC, type = "HC2")
welfare.H3.conservative.CI.partial <- coefci(welfare.H3.lm.partial, vcov. = vcovHC, type = "HC2")
pooled.H3.conservative.CI.partial <- conf_int(pooled.H3.lm.partial, 
                                              vcov = vcovCR(pooled.H3.lm.partial, 
                                                            cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                            type = "CR2"))

# effect of candidate gender for the liberal condition
DPRK.H3.lm.partial.reversed <- lm(DPRK.outcome ~ DPRK.female + I(1 - DPRK.liberal) + DPRK.aggressive + 
                                    order + DPRK.female : I(1 - DPRK.liberal), 
                                  data, subset = DPRK.correct == 1)
DPRK.H3.liberal.CI.partial <- coefci(DPRK.H3.lm.partial.reversed, vcov. = vcovHC, type = "HC2")
welfare.H3.lm.partial.reversed <- lm(welfare.outcome ~ welfare.female + I(1 - welfare.liberal) + welfare.aggressive + 
                                       order + welfare.female : I(1 - welfare.liberal), 
                                     data, subset = welfare.correct == 1)
welfare.H3.liberal.CI.partial <- coefci(welfare.H3.lm.partial.reversed, vcov. = vcovHC, type = "HC2")
pooled.H3.lm.partial.reversed <- lm(outcome ~ female + I(1 - liberal) + aggressive + order + 
                                      female : I(1 - liberal) + welfare, 
                                    pooled.data, subset = correct == 1)
pooled.H3.liberal.CI.partial <- conf_int(pooled.H3.lm.partial.reversed, 
                                         vcov = vcovCR(pooled.H3.lm.partial.reversed, 
                                                       cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                       type = "CR2"))

cairo_pdf("Figure_A4.pdf", width = 5.9, height = 3, pointsize = 8)
par(mar = c(5, 1, 1, 1))
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.8, 1.1), ylim = c(0, 14), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = 0, lty = 3, col = "gray")
points(DPRK.H1.lm.partial$coefficients[2], 13, pch = 19)
segments(DPRK.H1.CI.partial[2, 1], 13, DPRK.H1.CI.partial[2, 2], 13)
points(DPRK.H3.lm.partial.reversed$coefficients[2], 12, pch = 19)
segments(DPRK.H3.liberal.CI.partial[2, 1], 12, DPRK.H3.liberal.CI.partial[2, 2], 12)
points(DPRK.H3.lm.partial$coefficients[2], 11, pch = 19)
segments(DPRK.H3.conservative.CI.partial[2, 1], 11, DPRK.H3.conservative.CI.partial[2, 2], 11)
segments(DPRK.H3.liberal.CI.partial[2, 2] + 0.025, 12, DPRK.H3.liberal.CI.partial[2, 2] + 0.1, 12, lwd = 0.5)
segments(DPRK.H3.conservative.CI.partial[2, 2] + 0.025, 11, DPRK.H3.liberal.CI.partial[2, 2] + 0.1, 11, lwd = 0.5)
segments(DPRK.H3.liberal.CI.partial[2, 2] + 0.1, 12, DPRK.H3.liberal.CI.partial[2, 2] + 0.1, 11, lwd = 0.5)
text(DPRK.H3.liberal.CI.partial[2, 2] + 0.1, 11.5, 
     paste0(sprintf("%4.2f", DPRK.H3.lm.partial$coefficients[6]), " [\u2212", 
            sprintf("%4.2f", abs(DPRK.H3.conservative.CI.partial[6, 1])), ", ", 
            sprintf("%4.2f", DPRK.H3.conservative.CI.partial[6, 2]), "]"), cex = 0.8, pos = 4)
points(welfare.H1.lm.partial$coefficients[2], 8, pch = 19)
segments(welfare.H1.CI.partial[2, 1], 8, welfare.H1.CI.partial[2, 2], 8)
points(welfare.H3.lm.partial.reversed$coefficients[2], 7, pch = 19)
segments(welfare.H3.liberal.CI.partial[2, 1], 7, welfare.H3.liberal.CI.partial[2, 2], 7)
points(welfare.H3.lm.partial$coefficients[2], 6, pch = 19)
segments(welfare.H3.conservative.CI.partial[2, 1], 6, welfare.H3.conservative.CI.partial[2, 2], 6)
segments(welfare.H3.liberal.CI.partial[2, 2] + 0.025, 7, welfare.H3.conservative.CI.partial[2, 2] + 0.1, 7, lwd = 0.5)
segments(welfare.H3.conservative.CI.partial[2, 2] + 0.025, 6, welfare.H3.conservative.CI.partial[2, 2] + 0.1, 6, lwd = 0.5)
segments(welfare.H3.conservative.CI.partial[2, 2] + 0.1, 7, welfare.H3.conservative.CI.partial[2, 2] + 0.1, 6, lwd = 0.5)
text(welfare.H3.conservative.CI.partial[2, 2] + 0.1, 6.5, 
     paste0("\u2212", sprintf("%4.2f", abs(welfare.H3.lm.partial$coefficients[6])), " [\u2212", 
            sprintf("%4.2f", abs(welfare.H3.conservative.CI.partial[6, 1])), ", ", 
            sprintf("%4.2f", welfare.H3.conservative.CI.partial[6, 2]), "]"), cex = 0.8, pos = 4)
points(pooled.H1.lm.partial$coefficients[2], 3, pch = 19)
segments(pooled.H1.CI.partial[2, 4], 3, pooled.H1.CI.partial[2, 5], 3)
points(pooled.H3.lm.partial.reversed$coefficients[2], 2, pch = 19)
segments(pooled.H3.liberal.CI.partial[2, 4], 2, pooled.H3.liberal.CI.partial[2, 5], 2)
points(pooled.H3.lm.partial$coefficients[2], 1, pch = 19)
segments(pooled.H3.conservative.CI.partial[2, 4], 1, pooled.H3.conservative.CI.partial[2, 5], 1)
segments(pooled.H3.liberal.CI.partial[2, 5] + 0.025, 2, pooled.H3.conservative.CI.partial[2, 5] + 0.1, 2, lwd = 0.5)
segments(pooled.H3.conservative.CI.partial[2, 5] + 0.025, 1, pooled.H3.conservative.CI.partial[2, 5] + 0.1, 1, lwd = 0.5)
segments(pooled.H3.conservative.CI.partial[2, 5] + 0.1, 2, pooled.H3.conservative.CI.partial[2, 5] + 0.1, 1, lwd = 0.5)
text(pooled.H3.conservative.CI.partial[2, 5] + 0.1, 1.5, 
     paste0(sprintf("%4.2f", abs(pooled.H3.lm.partial$coefficients[7])), " [\u2212", 
            sprintf("%4.2f", abs(pooled.H3.conservative.CI.partial[7, 4])), ", ", 
            sprintf("%4.2f", pooled.H3.conservative.CI.partial[7, 5]), "]"), cex = 0.8, pos = 4)
text(rep(-0.8, 3), c(14, 9, 4), c("DPRK", "Public assistance", "Pooled"), pos = 4)
text(rep(-0.75, 3), c(13, 8, 3), rep("Female"), pos = 4)
text(rep(-0.75, 3), c(12, 7, 2), rep("Female | liberal"), pos = 4)
text(rep(-0.75, 3), c(11, 6, 1), rep("Female | conservative"), pos = 4)
axis(1, at = seq(-0.2, 0.6, 0.2))
mtext("Effects", side = 1, at = 0.2, line = 3)
dev.off()

### testing Hypothesis 4
## DPRK question
# linear regression
DPRK.H4.lm.partial <- lm(DPRK.outcome ~ DPRK.female + DPRK.liberal + DPRK.aggressive + order + DPRK.female : DPRK.aggressive, 
                         data, subset = DPRK.correct == 1)
# testing coefficients with the robust standard errors
DPRK.H4.result.partial <- coeftest(DPRK.H4.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(DPRK.H4.result.partial[, 1:2], 2), round(DPRK.H4.result.partial[, 4], 3))
# number of observations
length(DPRK.H4.lm.partial$residuals)

## welfare question
# linear regression
welfare.H4.lm.partial <- lm(welfare.outcome ~ welfare.female + welfare.liberal + welfare.aggressive + order + welfare.female : welfare.aggressive, 
                            data, subset = welfare.correct == 1)
# testing coefficients with the robust standard errors
welfare.H4.result.partial <- coeftest(welfare.H4.lm.partial, vcov. = vcovHC, type = "HC2")
# summary
cbind(round(welfare.H4.result.partial[, 1:2], 2), round(welfare.H4.result.partial[, 4], 3))
# number of observations
length(welfare.H4.lm.partial$residuals)

# conditional effect
welfare.H4.lm.partial.reversed <- lm(welfare.outcome ~ welfare.female + welfare.liberal + I(1 - welfare.aggressive) + 
                                       order + welfare.female : I(1 - welfare.aggressive), 
                                     data, subset = welfare.correct == 1)
welfare.H4.result.partial.reversed <- coeftest(welfare.H4.lm.partial.reversed, vcov. = vcovHC, type = "HC2")
cbind(round(welfare.H4.result.partial.reversed[, 1:2], 2), round(welfare.H4.result.partial.reversed[, 4], 3))

## pooled analysis
# linear regression
pooled.H4.lm.partial <- lm(outcome ~ female + liberal + aggressive + order + female : aggressive + welfare, 
                           pooled.data, subset = correct == 1)
# testing coefficients with the robust standard errors
pooled.H4.result.partial <- coef_test(pooled.H4.lm.partial, 
                                      vcov = vcovCR(pooled.H4.lm.partial, 
                                                    cluster = pooled.data[pooled.data$correct == 1, "id"], 
                                                    type = "CR2"))
# summary
cbind(round(pooled.H4.result.partial[, 1:2], 2), p = round(pooled.H4.result.partial$p_Satt, 3))
# number of observations
length(pooled.H4.lm.partial$residuals)